
    ##############################################################################################
    # R code to replicate
    # Figure 3 of Measuring Foreign Policy Positions of Members of the US Congress
    ##############################################################################################

    #setwd("C:\\0Jeong\\1_Projects\\Foreign_Policy\\0_Estimating_FP_Ideologies\\2015_PSRM\\Replication\\Replication_Material\\Supplemental_Material")

    Data <- read.table("Figure3_Data1.txt", header=T, sep="") # Roll Call Data
    Time <- read.table("Figure3_Data2.txt", header=T, sep="") # Time Index to be fed into MCMCdynmicIRTi1
    info2 <- read.csv("Figure3_Data3.csv", header=T)          # Legislator Information: ICPSR_ID, congress, party ID, and Name (Adapted from Voteview.com data)
    NOM_Dem.medians <- read.csv("Figure3_Data4.csv", header=T)$Dem.median  # Democratic median score of DW-NOMINATE dimension 1 (computed from Voteview.com data)
    NOM_Rep.medians <- read.csv("Figure3_Data4.csv", header=T)$Rep.median  # Republican median score of DW-NOMINATE dimension 1 (computed from Voteview.com data)


    library(MCMCpack)

    out <- MCMCdynamicIRT1d(Data, 
                            item.time.map=Time$time,
                            mcmc=20000, burnin=1000, thin=20,       
                            verbose=1000, store.item=FALSE,         # I set "store.item=TRUE" to produce outputs for Figure 6 in Supplemental Files. For a quick result, set it FALSE.          
                            theta.constraints=list(A9369A="+",     # Strom Thurmond
                                                   A14105A="+",    # Jesse Helms
                                                   A14920A="-",    # John Kerry 
                                                   A10808A="-"))   # Edward Kennedy

 
    out_ideal <- out[seq(from=1, to=nrow(out), length=min(1000, nrow(out))) , 1:3350]   # select ideal points only
    N <- ncol(out_ideal)

    # to check if converged, select 5 chains randomly.
    postscript("Supplemental_Figure15.eps", height=11, width=7.5, horizontal=F)
    n.size <- 5
    par(mfrow=c(n.size,2))
    for(i in sample(1:3350, n.size, replace=F)){
    plot(density(out_ideal[,i]), main="", xlab="")
    title(main=paste("Ideal Point of", colnames(out_ideal)[i]), font.main=1)
    plot(out_ideal[,i], type="l", xlab="iteration", ylab="")
    title(main="Traceplot", font.main=1)
    }
    dev.off()

    ### Geweke Diagnostics 
    # Note Geweke stat does not work well when the autocorrelation is high. 
    # If so, increase mcmc iterations and the thinning interval (say, from 10 to 50 or 100). (but it means longer simulations).
    # Also, make sure to change "store.item=TRUE" to FALSE to speed up simulation. 

    Geweke.stat1 <- rep(NA, N)
    for(i in 1:N){
    Geweke.stat1[i] <- geweke.diag(out_ideal[,i], frac1=0.4, frac2=0.6)$z
    }
    postscript("Supplemental_Figure16.eps", height=7, width=11, horizontal=F)
    hist(Geweke.stat1, breaks=100, xlim=c(-5, 5), xlab="Geweke Statistics", main="")
    abline(v=c(-1.96, 1.96), lty=2)
    dev.off()


